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Abstract 

The adjoint Fokker-Planck equation method is applied to study the runaway probability function 
and the expected slowing-down time for highly relativistic runaway electrons, including the loss of 
energy due to synchrotron radiation. In direct correspondence to Monte Carlo simulation methods, 
the runaway probability function has a smooth transition across the runaway separatrix, which can 
be attributed to effect of the pitch angle scattering term in the kinetic equation. However, for the 
same numerical accuracy, the adjoint method is more efficient than the Monte Carlo method. The 
expected slowing-down time gives a novel method to estimate the runaway current decay time in 
experiments. A new result from this work is that the decay rate of high energy electrons is very 
slow when E is close to the critical electric field. This effect contributes further to a hysteresis 
previously found in the runaway electron population. 


1 



It is well known that under an external electric field a superthermal electron in plasma 
can run away from the bounds of the collisional force and get accelerated to very high 
energy [T]. Runaway electrons can be produced in tokamak disruptions, which is an im¬ 
portant issue for disruption mitigation on the International Thermonuclear Experimental 
Reactor (ITER). Further studies have thus been motivated to address the dynamics of the 
runaway electrons in momentum space [2H1]. Recently it has been shown that in the presence 
of a magnetic field, the synchrotron radiation reaction force can play an important role in 
the runaway electron dynamics jU-E]. Combined with pitch angle scattering, the radiation 
force can produce additional stopping power [B]. The resulting effects include increase of the 
critical electric field Eq (which will be above the Connor-Hastie field Ec)|5], IH], and modifi¬ 
cation of the runaway electron growth and decay rate. In addition, the stopping power can 
help form an “attractor” in electron momentum space, which can lead to a bump-on-tail 
distribution [9]. Taking these into account, simulations [10] have produced results that quali¬ 
tatively agree with experiment [TT] for the electric field above which x-ray signals indicate a 
runaway population is sustained. It is believed that other effects like bremsstrahlung radia¬ 
tion and magnetic fluctuations can also influence the runaway growth and decay. To better 
understand these effects it is very important to develop theoretical tools that complement 
numerical simulations and can provide deeper physical insight into the phase-space structure 
of runaway electrons. 

In this paper we study the runaway electron dynamics in momentum space by solving the 
adjoint Fokker-Planck equation (EPF) in momentum space, which is a general method that 
offers significant conceptual and computational advantages. The homogeneous adjoint FPE 
was first introduced to study the first passage problem[l2]. It has been applied to calculate 
the neutron generation probability [T3|, the response function of the current drive[Tl|, and the 
runaway probability |15]. Here we demonstrate that the adjoint FPE can be used not only to 
study the probability function and its moments, as is often done, but also to calculate subtler 
and experimentally relevant quantities like the slowing-down time for existing superthermal 
electrons using the nonhomogeneous form of the FPE (see Appendix]]). This method takes 
into account all the terms in the kinetic equation, and improves upon the test particle 
method which ignores the diffusion term[T6l |T7|. In addition, the adjoint method is much 
more efficient than the Monte Carlo method since it can provide detailed information in all 
of momentum space by solving a single partial differential equation (PDF) once. It can be 
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more advantageous to study certain physical effects using the adjoint method, because the 
adjoint FPE has a direct relation to the standard Fokker-Planck equation. 

We demonstrate that the the runaway probability function shows a transition layer in 
momentum space, which agrees with the separatrix found by the test particle method. 
However, due to the effect of pitch angle scattering, the layer of finite width provides a 
smooth transition rather than a discontinuous transition represented by a step function (in 
momentum space). The expected slowing-down time we calculate characterizes the runaway 
electron beam decay time, which gives a new perspective to the study of runaway current 
decay in both the quiescent runaway electron (QRE)[TT] and the plateau[TH] regimes. The 
result shows that the electric field must be well below Eq for significant decay to occur. 

In the established model of runaway electron dynamics, when E is larger than the critical 
electric field and the radiation effect is weak, electrons initially in the high energy regime 
can continue to be accelerated and run away. On the other hand, electrons initially in the 
low energy regime will be decelerated and fall back into the Maxwellian population. Thus 
the destinations of electrons in the long time limit depend on their initial momentum. The 
radiation force can be an additional source of stopping power, but it can only dominate the 
electric force in the very high energy regime when E is much larger than the critical electric 
field. The kinetic equation for relativistic electrons can be written as[71 El ITS] . 




Z + 1 d _ 

2 ae ^ ^ di 



( 1 ) 


where p is the electron momentum (normalized to rriec), ^ is the cosine of the pitch angle. 


Z is the ion effective charge, E = E/Ec where E^ is the Connor-Hastie critical electric 
field Ec = rieC^lnA/ (dvre^meC^) and In A is the Coulomb logorithm, t = t/r where r is the 


relativistic electron collision time r = ruec/ (E^e), fr = TrI t and is the timescale for the 

synchrotron radiation energy loss Tr = ^'neQm\(? j (e'^H^). 

In the adjoint method, we define P{po,^o) as the runaway probability function, which 
means the probability for an electron that is initially at {po,^o) to eventually run away. As 


shown in Appendix | P satisfies the homogeneous adjoint equation of Eq. ([^, 

£[P]+C [P] + 5 [P] + [P] = 0, 


( 2 ) 
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where the four terms represent, respectively, the parallel electric held force, the drag force in 
collision operator, the pitch angle scattering, and the synchrotron radiation reaction force. 
The boundary conditions of P are set as P{p = Pmin,0 = 0, P{p = Pmax,0 = 1, where 
Pmin and Pmax are two boundaries in momentum space that are located far from the transition 
region. (The solution is checked to be insensitive to the boundary locations) 

We solve Eq. (|^ numerically using the hnite difference method, which is similar to the 
numerical method in Ref. UHl Figure shows P for E/E(. = 6, Z = 1 and = 100. The 
separatrix calculated using the test particle method in Ref. [I7| is also shown for reference. 
Note that the separatrix lies in the transition region of P {P between 0 and 1). However, 
we now have a smooth function that transtions from 0 to 1 rather than a step function. The 
width of the transition region depends on the amplitude of the pitch angle scattering term, 
which increases with Z. This transition region is not captured in the test particle method. 
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FIG. 1. The runaway probability for E = 6Ec, Z = 1 and fr = 100. 9 is the pitch angle. The red 
dashed line is the separatrix calculated using the test particle method in Ref. m- 

Figure 1^ shows the value of P as a function of p for ^ = 1 in the transition region. The 
result is benchmarked with a Monte-Carlo simulation result, which is achieved by sampling a 
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FIG. 2. The runaway probability function P for ^ = 1 (blue line), compared with the Monte-Carlo 
simulation result (red dashed line) and the separatrix calculated by the test particle method (black 
vertical line). The standard error of the Monte Carlo result is shown as error bars. 

large number of electrons that start at one initial position and follow the equation of motion 
that corresponds to the FPE Eq. (|^. We then count the electrons that hit the low and 
high energy boundaries after a certain time. The two results are close. Note that unlike the 
Monte-Carlo method which can take signihcant computer time, our method is fast and only 
requires solving the PDE once to obtain the probability function. IZDl 

P has an increasingly sharp transition from near zero to near unity as E increases to 
large values, which indicates that an electron with initial momentum above the separatrix 
is very likely to run away. However, as E decreases, the value of P above the separatrix 
reduces and eventually approaches zero. For E sufficiently small P becomes a flat function 
close to zero in most of momentum space, and only increases to 1 in a thin layer close to 
Pmax- fni This is because the electric held is so weak that it is dominated by the combination 
of the collisional drag and the radiation force. Thus, all electrons, regardless of their initial 
energy, will slow down to the low energy regime in a hnite time. This indicates the runaway 
population converts from growth to decay. 

In the decay phase the expected electron slowing-down time as a function of momentum, 
as opposed to the runaway probability, characterizes the decay of the runaway electron 
population. This can also be calculated using the adjoint method. Dehne T{p,^) as the 
expected time for an electron initially at (p, to reach the low energy boundary pmin or the 
high energy boundary Pmax- Note 1/T = l/Tg + 1/Tr, where Tg is the expected slowing down 
time and is the expected time to run away. The ratio of the two terms is (1 — P)/P. As 
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shown in Appendix [| T satisfies the nonhomogeneous adjoint FPE, 


8[T]+C[T]+S[T]+n[T] =-1. ( 7 ) 

The bonndary conditions are T{p = Pmin,0 = 0) = Pmax,0 = 0- Note that for 

E < Eq, the runaway probability is close to zero almost everywhere, so —)■ cxd and 

T ^ Ts, except for the region near the high energy boundary. 

Figure shows the calculated Ts{p,^) for ^ = 1 by solving Eq. 0, for E = 1 . 5 , Z = 1 . 
We see that Tg is a monotonically increasing function of p. For small radiation force (large 
fr) and E close to the critical field Eq, Tg has a large jump between the low and high energy 
regimes. This is because in the intermediate energy regime all the forces reach a balance and 
the motion is dominated by the diffusion effect, therefore electrons take a very long time to 
cross this barrier region through random walk. For large radiation force and E smaller than 
Eq this jump is very small or non-existent because the radiation force is strong and always 
dominates the electric field force. 

We also calculate the effect of energy loss due to large angle collisions by including a 
Boltzmann collision operator using the Mpller scattering cross sectionjH |22] in the adjoint 
equation. The result (dashed line in Figure]^ shows a significant decrease of slowing down 
time at the marginal case where E close to Eq, while for E much smaller than Eq the result 
does not change much. This decrease occurs because in the marginal case, the diffusion 
barrier formed by the balance of forces is significant. The large angle collisions, however, 
can help electrons cross the barrier since they can cause a high energy electron to lose a 
large fraction of its energy and fall directly into the low energy regime. However, the jump 
in Tg still exists. 

The expected slowing-down time can be used to estimate the runaway electron beam 
decay time in experiments, and help explain the runaway electron population hysteresis 
and distribution. In both the QRE and current quench experiments, due to the decreasing 
magnitude of E/E^, the runaway electron beam will have a transition from growth to decay. 
This means that at the beginning of the decay, there is already a population of high energy 
electrons formed by previous growth. The expected slowing-down time for these electrons 
determines the timescale for the runaway beam decay. In fact, if E is very close to Eq, 
the expected slowing-down time for the high energy electrons can be very long, due to the 
jump in Tg. This leads to a stagnation stage for the high energy electrons. The electric field 
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FIG. 3. The expected slowing-down time (normalized to r) as a function of p at .^ = 1 for 
E/Ec = 1.5, Z = 1 and 3 values of ?>. Dashed line is the results including large angle collision 
effects with In A = 30 

required for a significant decay to happen is thus far lower than the critical electric held 
an effect hrst captured by this model. This can contribute to a hysteresis [S] for the runaway 
electron population when the electric held is ramped up and down. Another indication from 
the Ts solution is that, due to very fast decay of the low energy electrons and extremely 
slow decay of the fast electrons, the electron population will tend to form a bump-on-tail 
distribution in the decay phase, which is diherent from the monotonic distribution in the 
growth phase[23]. Further investigations of the outcome of this distribution will be discussed 
in future work. 

Returning to the runaway probability P, as E is reduced P will suddenly change from a 
smooth transition across the separatrix to a hat function near zero, indicating the cessation 
of the runaway generation process. This sudden change in the structure of P can be used 
to determine the critical electric held Eq for runaway electron growth, which is above Ec in 
the presence of the radiation force. One should bear in mind that this critical value is not 
an absolute threshold, given that P always smoothly transitions to 1 at Pmax, in a thin layer 
near p^ax as E approaches Eq. We can, however, dehne a criterion based on the presence 
of a transition across the separatrix. For low Z the transition is rather abrupt, while for 
high Z it is smoother. Here we choose the a precise (but somewhat arbitrary) criterion that 
if P is above 0.005 in the region above the separatrix, which means an electron there has a 
0.5% probability to run to the high energy boundary, then the runaway generation process 
is active. We then obtain Eq from this criterion. 
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We have calculated Eq for 1 < Z < 10 and 10 < fr < 100, as shown in Figure The 
high energy boundary is chosen to be 30MeV. A convenient function that hts the result is 


X = 


Z + 1 


-^0 . V 

— = 1 + ax , 
Er 


a = 1.8587, 1 / = 0.6337. 


( 8 ) 


(9) 


(r^/r)3/4’ 

The error of the htting function is less than 5%. We have also benchmarked our result with 



FIG. 4. Eq (solid line) plotted as a function of A for various Z. Eq from Ref. [8] (dashed line) 
plotted for comparison. 


Eq in Ref. IH]in Figure]^ The two results are close for small Z, while for large Z our result is 
larger for small A but smaller for large fy. The discrepancy is mainly due to the difference 
in our dehnition of Eq and the uncertainty introduced by the smoother change of P for high 
Z. 

Note that this critical electric held may be different from experimental observations for 
several reasons. If the electron temperature is very low or the pitch angle scattering is strong, 
the critical energy required for an electron to run away (the least momentum above the 
separatrix) is very high, which results in a growth rate too low to be observed. Additionally, 
the observed electric held corresponding to the turning point of the signal in the QRE 
experiments [TT] can be higher than Eq, due to the energy dependence of the diagnostic and 
the redistribution of the runaway electron energy[7]. 

It is noteworthy that the the result of the adjoint method, especially the expected slowing- 
down time, depends highly on the energy dihusion mechanism in the kinetic model. In 
the results presented here, the model includes collisions and the synchrotron radiation re¬ 
action force. However, the method can be easily extended to study the influence of other 












physics on runaway electron dynamics, including Bremsstrahlung radiation[23] and magnetic 
fluctuations [23], by adding the corresponding operators into the adjoint FPE. In addition, 
the adjoint FPE can be used to study any dynamical system that has a separatrix or a sin¬ 
gular point, e.g. particle behavior close to the magnetic separatrix and the X-point. Future 
applications of this method to other areas are promising. 

We thank O. Embreus, I. Fernandez-Gomez, N. Fisch, T. Fiilop, P. Helander, E. Hirvijoki, 
J. Krommes, G. Papp and A. Stahl for useful discussions. The numerical calculations are 
conducted on the PPPL Beowulf cluster. This work is supported by the U.S. Department 
of Energy under Gontract No. De-FG02-03ER54696. 

Appendix: Adjoint Fokker-Planck equation. — Here we introduce the adjoint FPE 
method, begining with the homogeneous adjoint FPE. Gonsider a test particle in a sta¬ 
tionary stochastic system. Denote the particle’s coordinate as x with two boundaries Xmin 
and Xmax- The equation of motion of the test particle in the stochastic system can be de¬ 
scribed as 6x = x{t + 6t) —x{t) = v{x)5t + ^[x), where ^{x) is a random variable that satisfies 
(^(x)^(x)) = ^j2D{x)5t. The distribution function f{x,t) thus satisfies the FPE 



( 10 ) 


Define P(xo) as the probability of a test particle with initial coordinate x = Xq to first 
pass the boundary Xmax rather than Xmin- Note that because the system is stationary, P 
is time-independent. P(xo) can be expressed using the random walk probability density in 
terms of the particle’s next-step coordinate. 



G(xo, Sx)d6x, 


( 11 ) 


where G(xo, Sx) is the probability density that the coordinate can have a change of 6x in 
dt ii X = Xq, and we expand in powers of 6x in anticipation of taking the limit 5x —)■ 0. 
G(xo, Sx) satisfies the following properties 



( 12 ) 



( 13 ) 
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Using these equations, Eq. 0. and taking the limit 5t —)• 0, we obtain the differential 
equation for P(a;), 

r!P(P( 

(14) 


dPix) (PPix) 

n(x)—^ + D{x) , ’ = 0, 


dx 


dx"^ 


which is the adjoint equation of Eq. (10). According to the definition, the boundary con¬ 


ditions obeyed by P{x) are P{x = a^min) = 0, P{x = a^max) = 1- The probability for 
the particle to first pass the boundary at Xmin can be obtained simply from the relation 


Q{xq) = 1 — P(xo). Note that Q also satisfies Eq. (14). 

We next discuss the nonhomogeneous adjoint Fokker-Planck equation. Define T(a;o) as 
the expected time for a test particle that starts at x = Xq to first pass either boundary, Xmin 
or 3:^max- Similar to P, T can also be calculated through the random walk integral. 


T(xo) = J T{xo + 6x)G{xo, 5x)d6x + 6t 

= / [r(io) 


dT(xa) IdT(j'ci) ^ 


G(xo, 6x)d6x + 6t. 


Taking the limit St —)■ 0, P(x) is found to satisfy the differential equation 


(15) 


x(x) 


dT{x) 

dx 


D{x 


d^Tix) 


= - 1 , 


(16) 


which is analogous to Eq. (14) except the equation is now nonhomogeneous. The boundary 


conditions for T are simply T(x = Xmin) = 0, T(x = Xmax) = 0. 

Let us assume a particle source at x = xq. The rate for particles to pass one of the 
boundaries can be expressed as r = 1/T. Note that r = ri +r 2 , where ri and r 2 are the rate 
to pass the boundary at Xmin and x^ax- Both ri and r 2 can then be calculated according to 
the first passage probability, ri/r 2 = Q/P- 

It is important to point out that, though the derivation shown here is based simply on 
finding the adjoint FPE, the adjoint method can also be applied to more general kinetic 
PDFs. For example, to treat large angle collision effects one needs to use the Boltzmann 
collision operator, in which case the kinetic equation is an integro-differential equation rather 
than a differential one. However, one can still find the adjoint equation through integration 
by parts, or from the first line in Eq. ([IT| and Eq. (15) without performing Taylor expansion. 
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